/*************************************************************************
Purpose:
          Loads FES ascii data files, fills in NP hole by interpolation, 
          and writes it to disk using another name.
**************************************************************************    
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

int main(){  

  long i,j,l;
  model_s model;
  latlon_s rawfes;
  string basefolder(gaiafolder());
  
  /*
  string filename;
  double maxlat,minlat,maxlon,minlon;//Boundaries of actual model data.
  double deltalon, deltalat;//Lat/lon degrees in between model points.
  double mask;//Mask values occur on land.
  long long im,jm;//Number of data points in lon/lat directions.
  //int max_path=1024;//From FES loading code- max length of strings.
  //Used for loading FES amplitudes and phases along 2D grid.
  vector< vector<double> > amplitudes, phases;
  vector<double> lat,lon;//Defines coords of 2D grid for fes_amp/phase.
  //*/

  //Holds i,j indices of the hole at NP and the non-hole (edge) indices around it.
  vector<long long> hole_i,hole_j,edge_i,edge_j;
  model.options.basefolder = "data/FES2004/data/";
  model.options.type=1;
  basefolder.append(model.options.basefolder);

  define_model_titles_and_periods(model);
      
  for(l=0;l<(int)model.options.filenames.size();l++){
    hole_i.clear(); hole_j.clear();
    edge_i.clear(); edge_j.clear();

    rawfes.filename = basefolder;
    rawfes.filename.append(model.options.filenames.at(l));
    rawfes.discard_redundant_values = 1;
    read_rawfes(rawfes);
   
    cout<<"First, determine where the hole is, and what points lie on its edge."<<endl;
    //Begin search at end of latitude array, which is the north pole.
    int end = 0;
    j = (long long)rawfes.lat.size()-1;
    while(end==0 and j>0){//j cycles through latitude here.
      end=1;//Assume that hole is finished unless mask value is found at this lat.
      for(i=0;i<(long long)rawfes.lon.size();i++){//i cycles through longitude here.
        //If either amp or phase are the mask value, this is the hole.
        if(rawfes.amplitudes.at(j).at(i) == rawfes.mask or rawfes.phases.at(j).at(i) == rawfes.mask){
          hole_i.push_back(i); hole_j.push_back(j);
          end=0;//This is a hole value, so hole still exists at this lat.
        }
        else{
          edge_i.push_back(i); edge_j.push_back(j);
        }
      }
      j--;//Move south to the next latitude.
    }
    cout<<"Found "<<hole_i.size()<<" points in hole, and "<<edge_i.size()<<" points on the edge, with the southern-most hole point at latitude "<<rawfes.lat.at(hole_j.at(hole_j.size()-1))<<endl;
    
    cout<<"Filling hole by interpolating edge points."<<endl;
    //For every point in the hole, cycle through edge points and calculate
    //spherical-angle-weighted interpolated average using sin/cos rather than amp/phase.
    for(j=0;j<(long long)hole_j.size();j++){//here, j cycles through hole pts.
      double avg_sin = 0, avg_cos = 0, weight = 0;
      for(i=0;i<(long long)edge_i.size();i++){//i cycles through edge pts.
        double angle = spherical_distance(rawfes.lat.at(hole_j.at(j)),rawfes.lon.at(hole_i.at(j)),rawfes.lat.at(edge_j.at(i)),rawfes.lon.at(edge_i.at(i)));
        weight += 1/angle;//Record total weights for normalization later.
        avg_sin += (rawfes.amplitudes.at(edge_j.at(i)).at(edge_i.at(i))*sin(rawfes.phases.at(edge_j.at(i)).at(edge_i.at(i))*deg2rad))/angle;//Weight according to inverse separation angle. 
        avg_cos += (rawfes.amplitudes.at(edge_j.at(i)).at(edge_i.at(i))*cos(rawfes.phases.at(edge_j.at(i)).at(edge_i.at(i))*deg2rad))/angle;//Weight according to inverse separation angle. 
      }
      avg_sin /= weight; avg_cos /= weight;
      rawfes.amplitudes.at(hole_j.at(j)).at(hole_i.at(j)) = sqrt(pow(avg_sin,2)+pow(avg_cos,2));
      rawfes.phases.at(hole_j.at(j)).at(hole_i.at(j)) = rad2deg*atan2(avg_sin,avg_cos);
    }

    rawfes.filename = basefolder;
    rawfes.filename.append("new/");
    rawfes.filename.append(model.options.filenames.at(l));
    write_rawfes(rawfes);
  }
  
  return 0;
}
